#include <sstream>

TGraph2D* sigplot_2D(TString plot, int gggg, int tau=2,
		     TString text="Pythia",
		     TString xaxis="#tau_{2}/#tau_{1}", 
		     TString yaxis="#Sigma_{J} m_{J}",
		 double cxn_QCD=7.60016e9, 
		 double cxn_glu1=0.287e8, 
		 double cxn_glu2=92.6e3, 
		 double cxn_glu3=1.01e3, 
		 double cxn_glu4=.435e2)
{
  TH2* psig=(TH2*) _file1->Get(plot);
  TH2* pbkg=(TH2*) _file0->Get(plot);
  TH2* pbkg_stat = (TH2*) _file2->Get(plot);

  psig->RebinX(5);
  psig->RebinY(5);
  pbkg->RebinX(5);
  pbkg->RebinY(5);

  pbkg_stat->RebinX(5);
  pbkg_stat->RebinY(5);


  if(tau==3)
    xaxis="#tau_{3}/#tau_{1}";

  if(tau==32)
    xaxis="#tau_{3}/#tau_{2}";
 
  float nx=psig->GetNbinsX();
  float ny=psig->GetNbinsY();

  float increx =psig->GetXaxis()->GetXmax()/nx;
  float increy =psig->GetYaxis()->GetXmax()/ny;

  TH1D* qcd_nevt= (TH1D*) _file0->Get("nevt");
  TH1D* glu_nevt= (TH1D*) _file1->Get("nevt");


  float cxn_glu=cxn_glu1;

  if(gggg == 300)
    cxn_glu=cxn_glu2;
  if(gggg == 600)
    cxn_glu=cxn_glu3;
  if(gggg == 900)
    cxn_glu=cxn_glu4;

  psig->Scale(cxn_glu / glu_nevt->Integral(0,1) );
  pbkg->Scale(cxn_QCD / qcd_nevt->Integral(0,1) );

  //psig->Scale(lumi*cxn/nevt);

  float cutx=increx;
  
  vector<double> v_cutx, v_cuty, v_sig, v_stat;
  double opt_cutx=1.0;
  double opt_cuty=0;
  double opt_sig=0;
  double opt_stat=0;
  double opt_nsig=0;
  double opt_nbkg=0;

  for(int i=0; i<nx; i++)
    {
      float cuty=increy;
      for(int j=0; j<ny; j++)
	{
	  double c_sig=psig->Integral(0, i+1, j+2, ny+1);
	  double c_bkg=pbkg->Integral(0, i+1, j+2, ny+1);
	  double c_bkg_stat=pbkg_stat->Integral(0, i+1, j+2, ny+1);

	  double c_num=0;

	  if(c_bkg_stat <10)
	    c_num=0;
	  else
	    c_num=(c_sig/sqrt(c_bkg));

	  cout<<"Cut: "<<cutx<<","<<cuty;
	  cout<<"\tSig: "<<c_num<<","<<c_sig<<","<<c_bkg;
	  cout<<"\tStat: "<<c_bkg_stat<<" evts"<<endl;

	  if(c_bkg<10)
	    continue;
	  
	  v_sig.push_back(c_num);
	  v_cutx.push_back(cutx);
	  v_cuty.push_back(cuty);
	  
	  v_stat.push_back(c_bkg_stat);

	  //cout<<cutx<<","<<cuty<<","<<c_sig/sqrt(c_bkg)<<endl;
	  cuty+=increy;
	  if(c_num > opt_sig)
	    {
	      opt_nsig=c_sig;
	      opt_nbkg=c_bkg;
	      opt_cutx=cutx;
	      opt_cuty=cuty;
	      opt_stat=c_bkg_stat;
	      opt_sig=c_num;
	    }

	}
      cutx+=increx;
    }

  cout<<"Opt: "<<opt_cutx<<","<<opt_cuty;
  cout<<"\tSig: "<<opt_sig<<","<<opt_nsig<<","<<opt_nbkg;
  cout<<"\tStat: "<<opt_stat<<" evts"<<endl;

  stringstream sout;
  sout<<text<<" "<<gggg<<" Gluino S/#sqrt{B} vs. "<<xaxis
      <<" and "<<yaxis<<" cuts";

  stringstream sout1;
  sout1<<text<<" "<<gggg<<" Gluino unweight number of events vs. "<<xaxis
      <<" and "<<yaxis<<" cuts";

  TGraph2D* a=new TGraph2D(v_cutx.size(), &v_cutx[0], &v_cuty[0], &v_sig[0]);
  a->SetTitle(sout.str().c_str());

  //TGraph2D* b=new TGraph2D(v_cutx.size(), &v_cutx[0], &v_cuty[0], &v_stat[0]);
  //b->SetTitle("QCD number of events");

  if(tau==2)
    a->SetName("tau21");
  if(tau==3)
    a->SetName("tau31");
  if(tau==32)
    a->SetName("tau32");
  
  
  a->GetXaxis()->SetTitle(xaxis);
  a->GetXaxis()->SetTitleOffset(1.0);
  a->GetYaxis()->SetTitleOffset(1.5);
  a->GetYaxis()->SetTitle(yaxis);
  a->GetZaxis()->SetNdivisions(505, kTRUE);

  /*
  b->GetXaxis()->SetTitle(xaxis);
  b->GetXaxis()->SetTitleOffset(2);
  b->GetYaxis()->SetTitleOffset(2);
  b->GetYaxis()->SetTitle(yaxis);
  */


  TCanvas* c1=new TCanvas("c1", "c1", 600, 600);
  a->Draw("colz");


  TString save=text;


  if(tau==2)
    {
      save+="_";
      save+=gggg;
      save+="_1D_tau21.png";
    }

  if(tau==3)
    {
      save+="_";
      save+=gggg;
      save+="_1D_tau31.png";
    }
  
  if(tau==32)
    {
      save+="_";
      save+=gggg;
      save+="_1D_tau32.png";
    }
  

  c1->SaveAs(save);

  //  TCanvas* c2=new TCanvas("c2", "c2", 800, 600);
  //  b->Draw("TRI1");
  
  return a;
}
